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ABSTRACT 

The MCPep server (http://bental.tau.ac.il/MCPep/) is 
designed for non-experts wishing to perform Monte 
Carlo (MC) simulations of helical peptides in associ- 
ation with lipid membranes. MCPep is a web imple- 
mentation of a previously developed MC simulation 
model. The model has been tested on a variety of 
peptides and protein fragments. The simulations 
successfully reproduced available empirical data 
and provided new molecular insights, such as the 
preferred locations of peptides in the membrane 
and the contribution of individual amino acids to 
membrane association. MCPep simulates the 
peptide in the aqueous phase and membrane envir- 
onments, both described implicitly. In the former, 
the peptide is subjected solely to internal conform- 
ational changes, and in the latter, each MC cycle 
includes additional external rigid body rotational 
and translational motions to allow the peptide to 
change its location in the membrane. The server 
can explore the interaction of helical peptides of 
any amino-acid composition with membranes of 
various lipid compositions. Given the peptide's 
sequence or structure and the natural width and 
surface charge of the membrane, MCPep reports 
the main determinants of peptide-membrane inter- 
actions, e.g. average location and orientation in the 
membrane, free energy of membrane association 
and the peptide's helical content. Snapshots of 
example simulations are also provided. 

INTRODUCTION 

Biological membranes separate between the interior 
of cells and organelles and the exterior environment. 



The membranes are composed mostly of lipid and 
protein molecules of various types (1). Roughly 
speaking, the lipid molecules are packed heterogeneously 
and organized in a fluid mosaic, mutually affecting each 
other. Biological membranes and protein- and peptide- 
lipid interactions have been investigated extensively with 
a wide range of biophysical and computational techniques 
(2). The goal of these studies has been to decipher the 
mechanism of action of membrane-active peptides, such 
as antimicrobial and viral-fusion peptides, and to deduce 
general principles of folding, energetics, stability, structure 
and function of membrane proteins. 

Several computational approaches are commonly 
applied for investigation of protein- and peptide-lipid 
interaction (3-5). Our group developed a model that 
performs Monte Carlo (MC) simulations of the inter- 
action of helical peptides with lipid bilayers. The simula- 
tion model has been tested on a variety of antimicrobial 
peptides, including magainin2 (6), penetratin (6), melittin 
(7), NKCS (8) and novicidin (9). These studies have fur- 
thered our understanding of the action mechanisms and 
selectivity of the antimicrobial peptides. The simulations 
correlated very well with the empirical data, and in some 
cases, they were used to guide further experimental effort 
(7-9). Additionally, the model was used to investigate the 
channel-forming peptide M28 from the acetylcholine 
receptor 5-subunit (10). The energetically preferable con- 
figuration and the 3-dimensional (3D) structure of the 
M2§ peptide in the membrane were correctly predicted. 
We also used the model to predict stable conformations 
of the outer membrane 45 (OM45) protein on interaction 
with the mitochondrial outer membrane (11). 

Here, we present MCPep (http://bental.tau.ac.il/ 
MCPep/): a web implementation of the MC simulation 
model. MCPep can be used to study short membrane- 
active peptides, such as antimicrobial and fusion 
peptides, short membrane proteins (see the Case Study 
described later in the text) or fragments of large 
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membrane proteins. The server takes as input the peptide's 
sequence or structure, fractions of acidic and zwitterionic 
lipids in the membrane, the native width of the hydrocar- 
bon region of the membrane and the ionic strength of the 
aqueous solution. The output includes the peptide helicity 
and average orientation in the membrane, the free energy 
of its membrane-association and the decomposition of the 
free energy. MCPep also presents snapshots of example 
simulations. In addition, the conformations obtained 
from the simulations are clustered, and the centroid con- 
formations of the largest clusters are provided in Protein 
Data Bank (PDB) format. 

THE MC MODEL 

The peptide is described in a reduced way; each amino 
acid is represented by two interaction sites, corresponding 
to its a-carbon and side chain (10). The interaction sites of 
sequential a-carbons are connected by virtual bonds. The 
membrane hydrophobicity is represented as a smooth 
profile, corresponding to the native thickness of the hydro- 
carbon region (10). 

The total free energy difference between a peptide in the 
aqueous phase and in the membrane (/1G tota i) is calculated 
as (12,13) follows: 

AGtotal = AGcon + AGdef+AGcoul+AGsol+AGimm + AGlip 

(1) 

where ^G con is the free energy change owing to 
membrane-induced conformational changes in the 
peptide. At constant temperature (T), it can be calculated 
as AGcon = AE — TAS, where AE is the internal energy 
difference between the peptide in the aqueous phase and in 
the bilayer. The internal energy is a statistical potential 
derived from available 3-dimensional (3D) protein struc- 
tures (14). The energy function assigns a score (energy) to 
each peptide conformation according to its abundance in 
the PDB. Common conformations are assigned high 
scores (low energy), whereas rare conformations are 
assigned lower scores (higher energy). zlS refers to the 
entropy difference between the aqueous solution and 
membrane-bound states, whereas the entropy (S) in each 
state is determined by the distribution of the virtual bond 
rotations in the reduced peptide representation. 

zlGdef is the free energy penalty associated with fluctu- 
ations of the membrane thickness around its native 
(resting) value, calculated following the estimation of 
Fattal and Ben-Shaul (15). Their calculations were based 
on a statistical-thermodynamic molecular model of the 
lipid chains. Their model fits a harmonic potential of the 
form: AGdef= cdAL 2 , where Ah is the difference between 
the native and actual width of the membrane; m is a 
harmonic-force constant related to 9 the membrane 

elasticity and is equal to co = 0.22kT/A (15). 

zlGcoui stands for the Coulombic interactions between 
titratable residues of the peptide and the (negative) 
surface charge of the membrane. It is calculated 
using the Gouy-Chapman theory that describes how the 
electrostatic potential depends on the distance from the 
membrane surface in an electrolyte solution (6). 



AG so \ is the free energy of transfer of the peptide from 
the aqueous phase into the membrane. It takes into 
account electrostatic contributions resulting from 
changes in solvent polarity, as well as nonpolar effects, 
both resulting from differences in van der Waals inter- 
actions of the peptide with the membrane and aqueous 
phases and from solvent structure effects. z1Gj mm is the 
free energy penalty resulting from the confinement of the 
external translational and rotational motion of the peptide 
inside the membrane. AGy ip is the free energy penalty re- 
sulting from the interference of the peptide with the con- 
formational freedom of the aliphatic chains of the lipids 
in the bilayer while the membrane retains its native 
width. The latter three terms, i.e. z1G so i, AG[ mm , AGn p , 
collectively referred to as z!G S il, are calculated using the 
Kessel and Ben-Tal hydrophobicity scale (12). The scale 
accounts for the free energy of transfer of the amino acids, 
located in the center of a polyalanine a-helix, from the 
aqueous phase into the membrane midplane. To avoid 
the excessive penalty associated with the transfer of 
charged residues into the bilayer, in the model, the titrat- 
able residues are neutralized gradually on insertion into 
the membrane, so that a nearly neutral form is desolvated 
into the hydrocarbon region. 

To calculate the free energy change in Equation 1, the 
peptide is simulated in the aqueous phase and membrane 
environment. In the aqueous phase, the peptide is sub- 
jected only to internal conformational modifications. In 
one MC cycle, the number of internal modifications per- 
formed is equal to the number of residues in the peptide. 
In the membrane, each MC cycle includes additional 
external rigid body rotational and translational motions 
to allow the peptide to change its location in the 
membrane and its orientation with respect to the 
membrane normal. A standard MC protocol is used, 
and acceptance of each move is based on the Metropolis 
criterion and the free energy difference between the new 
and old states (16). 



THE WEB SERVER 

Input 

The user is asked to provide the peptide sequence or 
upload an initial 3D structure in PDB format. The user 
should also specify the fraction of acidic lipids (mol 
percent) in the membrane and ionic strength of the 
aqueous solution. Additionally, the user may set the 
native width of the hydrocarbon core of the membrane 
according to the lipid types. To assist the user in doing 
so, the web server provides a link to recent X-ray studies 
of phospholipids [reviewed in reference (17)]. The user can 
also specify an e-mail address where a link to the results 
page will be sent once the simulations are completed 
(optional). 

Processing 

If a peptide sequence is provided as initial input, an initial 
canonical a-helix (phi = —58, psi = —47) model structure 
of the peptide is constructed, with side chains modeled 
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Table 1. Thermodynamic characteristics for syb2(75-l 16) in TM and 
surface configurations 
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AG totaI (kT) 


-53.2 ± 0.3 


-44.4 ± 0.3 


AG Con (kT) 


-1.5 ± 0.3 


5.3 ± 0.3 


TAS (kT) 


-32.0 ± 0.0 


-33.1 ± 0.3 


AE; nl (kT) 


-33.5 ± 0.3 


-27.8 ± 0.3 


AGsil (kT) 


-40.8 ± 0.1 


-37.9 ± 0.1 


AG def (kT) 


0.4 ± 0.00 


0.3 ± 0.0 


AGcoui (kT) 


-11.3 ± 0.0 


-12.1 ± 0.0 


M widlh (A) 


27.6 ± 0.0 


30.1 ± 0.0 


^center Ov 


10.6 ± 0.0 


18.7 ± 0.0 


Tilt (°) 


29.4 ± 0.2 


73.2 ± 0.2 



All values are reported as means ± standard error. The free energy 
terms are defined in Equation 1. M widlh : the width of the membrane 
hydrophobic core. Z cenler : the average distance of the peptide's center of 
mass from the membrane midplane. The Z-axis is the membrane 
normal, and the origin coincides with the membrane midplane. Tilt: 
the angle between the N'-to-C vector of the peptide's helical core 
and membrane normal. 

using Scwrl 4.0 (18); otherwise, the structure provided by 
the user is taken as the initial conformation. For simula- 
tions in the aqueous phase (i.e. without the membrane), 
the initial structure is used as it is. The simulations 
are carried out in three to five independent runs of 
500000-900 000 MC cycles each; the recommended 
number and length of the runs are set as defaults but the 
user can use alternatives. In a membrane, a helical peptide 
typically adheres to one of the two following configur- 
ations: transmembrane (TM) orientation, with the princi- 
pal axis of the helix positioned approximately along the 
membrane normal; or surface orientation, with the axis 
approximately parallel to the membrane surface. The tran- 
sition between the two configurations is associated with a 
high desolvation free energy barrier (2). Thus, for simula- 
tions in the membrane environment, each of the two con- 
figurations is used as the initial orientation. The number of 
independent simulations in each configuration and the 
number of cycles in each individual simulation are the 
same as in simulations in the aqueous phase. However, 
it should be noted that a simulation that starts in TM 
orientation may end in surface orientation if the peptide 
is amphipathic, too short to span the bilayer or too polar. 
The transition from a surface to TM orientation is highly 
unlikely owing to the high free energy penalty associated 
with the insertion of the polar backbone of the helix 
termini into the membrane. 

Output 

The web server's output includes the helical content of the 
peptide, i.e. the total number of residues that are in helix 
conformation, in the aqueous phase and in the membrane, 
calculated as described in reference (19). A residue is 
defined as being in a helical conformation if the rotational 
angle of each of its two adjacent bonds lies within the 
interval typical for a helix (-120° ± 30°) (14). As the ro- 
tational angles of the peptide's termini cannot be defined, 
the two residues at the N-terminus and the two residues at 
the C-terminus are neglected. Energetically favorable 
orientations of the peptide in the membrane are 



reported, including the average locations of the peptide's 
a-carbons. The free energy of membrane association and 
the decomposition of the free energy (Equation 1) are also 
presented. The membrane thickness is allowed to vary 
from its native value by up to 20%. The (perturbed) thick- 
ness along the simulation is also provided. In addition, 
snapshots of example simulations, both in the aqueous 
phase and in the membrane, are presented. The conform- 
ations obtained from the simulations are clustered with a 
preset cut-off for the root mean square deviation between 
the a-carbons. The default cut-off value is 3 A, but the 
user is free to alter the value. The centroid conformation 
of the largest cluster from each simulation is reported 
(in PDB format), along with its average orientation in 
the membrane. 



CASE STUDY 

The synaptic vesicle protein synaptobrevin2 (syb2, also 
known as vesicle-associated membrane protein 2 
(VAMP2)) is part of the so-called Soluble N- 
ethylmaleimide sensitive fusion protein Attachment 
Protein Receptor (SNARE) complex, responsible for 
fusion of neuronal vesicles with the presynaptic 
membrane in the neuron. Although the process is essential 
for synaptic neurotransmitter release, the exact mechan- 
ism of vesicle fusion is not yet clear. Here, we studied the 
membrane interaction of syb2 to illustrate the 
functionalities of the MCPep server. We used a peptide 
corresponding to the sequence of the putative TM and 
juxtamembrane regions of the protein, namely, residues 
75-116 (UniProt id P63027). The fraction of the charged 
lipids in the membrane was set to 30%, similar to that of 
neuronal vesicles (20). The results showed that both TM 
and surface orientations of syb2 in the membrane are 
stable, with the former being somewhat more favorable 
energetically (Table 1). The helical content of the 
peptide increased on interaction with the membrane, as 
anticipated (Figure 1). The average TM and surface orien- 
tations of syb2, as predicted by MCPep, are presented in 
Figure 2. In both surface and TM configurations, the 
hydrophobic residues are inserted into the membrane 
core, whereas the hydrophilic residues partition into the 
aqueous phase (Figures 2 and 3). 

We suggest that the exceptionally high stability of the 
TM configuration (— 53kT; Table 1) is required for syb2's 
function. Syb2 is located in vesicles, whereas its SNARE- 
complex partners, i.e. syntaxin and Soluble N- 
ethylmaleimide sensitive fusion protein Attachment 
Protein (SNAP)-25, are located in the presynaptic 
plasma membrane (21,22). Syntaxin and two molecules 
of SNAP-25 associate with syb2 via a process known as 
zippering, which starts at the mostly unstructured 
N-termini of the proteins. On association, the proteins 
form a remarkably stable tight bundle of four parallel 
coiled-coil a-helices, where each helix is from a different 
SNARE partner. As zippering continues towards the 
C-termini, syb2 pulls the vesicle towards the plasma 
membrane (21,22). This process eventually leads to the 
fusion of vesicular and plasma lipid bilayers. Clearly, 
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Figure 1. The average helical content of syb2 versus residue number as 
predicted by the MCPep server, (a) syb2 in the aqueous phase, (b) syb2 
in TM configuration within a membrane containing 30% anionic lipids. 
The helical content increases on interaction with the membrane, as 
anticipated. 
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Figure 2. The average location of the amino acids of syb2 in the 
membrane in (a) surface and (b) TM orientations. The membrane 
includes 30%-charged lipids. The horizontal dashed lines designate 
the location of the phosphate groups of the lipid polar heads. The 
hydrophobic residues (G, A, V, L, I, F) are in orange, polar residues 
(M, C, T, S, W, Y, H, Q, N) in green and charged amino acids 
(R, K, E) in blue. 
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Figure 3. The centroid conformation of the largest cluster of syb2 in 
TM orientation. The peptide is colored according to the hydrophobicity 
scale. The average location of the phosphate heads of the two leaflets of 
the lipid bilayer is represented by the red rectangles. 



vesicle pulling can only occur if both syb2-vesicle associ- 
ation and syb2/syntaxin/SNAP-25 interaction are strong. 
Indeed, the stability free energy of the SNARE complex is 
estimated to be between 18 (23) and 35 (24) kT, which is 
indicative of the strong forces required for membrane 
fusion to occur (21). In this respect, the strong anchoring 
of syb2 to the membrane (— 53kT; Table 1) is compatible 
with these estimates. For comparison, we also simulated 
the interaction of the putative TM region of syntaxin 
(UniProt id Q16623, residues 247-288) with a membrane 
containing 30% charged lipids and obtained a similar 
value of — 47kT. The free energy of membrane association 
of both proteins is comparable with the free energy 
required for hemi-fusion, that is the fusion of the outer 
leaflets of the bilayers, estimated at about 40-50kT (24). 



DISCUSSION 

MCPep allows rapid simulations of helical peptides in as- 
sociation with biological membranes. To minimize the 
computational burden, the server relies on reduced repre- 
sentation of the peptide. For similar reasons, it uses an 
implicit model of the membrane, which is described as a 
hydrophobic profile of preset thickness, and a surface 
charge (6,10). Thus, the computational model captures 
only certain characteristics of the complex peptide-lipid 
system. Clearly, other properties that might affect 
peptide-lipid interactions are missing in the model. One 
of these is the membrane curvature. The flat representa- 
tion of the membrane surface can roughly describe a cell 
membrane or a large vesicle. However, the curvature of 
small vesicles or micelles, which are comparable in size to 
the peptide length, cannot be overlooked and might play a 
significant role in the interaction with peptides. Another 
example is the phase of the lipids, which affects the fluidity 
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and rigidity of the bilayer. The MC model describes 
lipids in their crystalline-liquid phase, which are more sus- 
ceptible to stretching and bending than gel-phase lipids 
(25). The liquid phase bilayer adjusts more easily to the 
presence of membrane proteins. Lipid de-mixing is 
another characteristic that is not incorporated into the 
model. In this process, which occurs in interactions with 
cationic peptides, negatively charged lipids migrate to the 
peptide-membrane interaction zone (26-28). As a conse- 
quence, the local charge density in the interaction zone 
increases, which may induce an increased local concentra- 
tion of the peptide on the membrane's surface. 

The reduced representation also precludes investigation 
of specific interactions. For instance, it does not explicitly 
take into account hydrogen bonds and salt bridges 
between the side chains of polar amino acids or between 
the amino acids and the lipids 1 polar heads. Thus, the 
model is not suitable to describe peptide-peptide inter- 
actions, and it does not explicitly account for packing of 
the lipid chains against the peptide. Yet, the use of MC 
(rather than molecular dynamics) simulations facilitates 
efficient sampling in conformational/configurational 
space, making this approach suitable for studying thermo- 
dynamic quantities. Indeed, the previous studies using our 
MC model have demonstrated that computations using 
the approximations described earlier in the text are 
feasible and that many important characteristics of 
peptide-membrane interactions are captured (6-11). 

Many of the observations derived from the simulations 
depend on the water-to-membrane partition free energy 
of the amino acids, i.e. the hydrophobicity scale used 
in MCPep. We used the Kessel and Ben-Tal scale, derived 
from continuum (implicit)-solvent model calculations 
of the water-to-membrane transfer free energy of 
polyalanine-based peptides (12). The calculations were 
based on PARameter SEt (PARSE), a parameter set that 
was derived specifically to reproduce the experimental par- 
tition free energy of small molecules between water and 
liquid alkanes (29). The continuum-solvent model with 
PARSE has been used successfully to study thermodynamic 
and kinetic aspects of the transfer of several peptides 
and protein fragments into the lipid bilayer (10,13,30-33). 
The Kessel and Ben-Tal hydrophobicity scale was subse- 
quently exploited within the MC simulation model (6-10). 
The scale's hydrophobicity values are significantly lower 
in magnitude than those of other scales because the Kessel 
and Ben-Tal scale includes the free energy penalty owing 
to the transfer of the polar helix backbone from the 
aqueous phase into the membrane core (12). The scale is 
also unique in its extreme asymmetry: the free energy 
penalty associated with the water-to-membrane transfer of 
a single highly polar residue, such as Lysine in its neutral 
form, is approximately three times larger in magnitude than 
the gain owing to the transfer of a highly hydrophobic 
residue, such as Leucine (12). 



PLANS 

The methodology described herein can potentially be 
applied towards predicting TM protein structures. The 



most accurate models of protein structures are achieved 
through homology modeling, in which a model is inferred 
from a known high-resolution structure of a homologous 
protein (34). A crucial step in the modeling process is 
alignment of the query and template sequences (35), 
required to identify the TM helices of the query protein 
and match them to the template. Alignment is especially 
challenging if the query and target proteins share low 
sequence identity (36). The MCPep server could be used 
for identification of TM helices. In addition, it could 
predict the membrane boundaries for existing structures 
of TM proteins. Further investigation is required to 
explore these possibilities. 



CONCLUSIONS 

The MCPep server is capable of exploring the interactions 
of helical peptides of up to 50 residues with membranes 
of various types. Its results provide valuable mechanistic 
insights into membrane-peptide interactions, as demon- 
strated by the case study. It is preferable, when possible, 
to correlate the results with empirical data, and, thus, the 
conditions simulated in calculations should be similar to 
those used in the experiments. Model parameters that can 
be matched to experimental conditions include the 
fraction of acidic lipids in the membrane, the ionic 
strength of the surrounding aqueous solution and the 
native width of the membrane hydrocarbon core. The 
server can be easily handled even by non-expert users. 
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